Priority Effects
K. Carter Wynne
Overview
The role of temporal dispersal patterns in building diverse tallgrass prairie plant communities
Background:
Establishing diverse plant communities is critical since diversity is linked to ecosystem health. However, recreating tallgrass prairies with high plant diversity has been challenging, especially for early-season species. Variation in species arrival may influence plant community composition and diversity. Since prairies exhibit temporal patterns of seed dispersal, planting all species simultaneously could forgo phenological differences that promote species coexistence. Therefore, we investigated whether manipulating plant species’ arrival according to natural dispersal phenology influences reconstruction outcomes.
In 2021, we manipulated the arrival of 36 species via seed additions of (i) species in the order of peak dispersal timing, (ii) summer dispersing species (first peak in dispersal activity before September 1st) followed by 18 fall dispersing species, (iii) 18 fall dispersing species followed by 18 summer dispersing species, (iv) all species simultaneously. Additionally, we had a negative control that had no seed additions. One year later, we found that differences in seeding treatment influenced the diversity and composition of reconstructed communities. Species arriving later had less cover than when seeded with priority, particularly for summer dispersing species. Overall, our study provides evidence of priority effects in reconstructed grasslands and suggests that the timing of seed additions post-disturbance influences restoration outcomes.
Author: Katherine Carter Wynne (Wynnekat@msu.edu)
Co-authors: Lauren L. Sullivan
Created: 06 December 2022
Files:
PE_Cover_Inner_Fall_2021_Cleaned.xlsx - vegetative cover taken at peak biomass in 2021 (August 27th - August 30th); only half of the lumped seedings occured.
PE_Cover_Inner_Summer_2022_Cleaned.xlsx - vegetative cover measured in Summer 2022 (June 27th - June 28th)
PE_Cover_Inner_Fall_2022_Cleaned.xlsx - vegetative cover taken at peak biomass in 2022 (August 25th - August 26th)
PE_Cover_Inner_Summer_2023_Cleaned.xlsx - vegetative cover measured in Summer 2023 (June 21st - June 23rd)
PE_Cover_Inner_Fall_2023_Cleaned.xlsx - vegetative cover taken at peak biomass in 2023 (August 10th - August 11th)
Bradford_10_Year_Weather.xlsx - weather (temperature and precipitation) data from the last 10 years at Bradford Research Farm
Priority_Effects_Fall_2022_Light.xlsx - light data from August 2022; measured at the same time as peak biomass cover
Priority_Effects_Fall_2023_Light.xlsx - light data from August 2023; measured at the same time as peak biomass cover
PE_Bare_Ground_Cleaned_2021.xlsx - Litter and bare ground cover from August 2021
Species_List.xlsx - Species list
R Version: R 4.2.2
RStudio Version: 2023.06.1+524
Package Version: Found in the .readme document in GitHub Repository
Last updated: 29 September 2023
Overall Setup
Materials and Methods
Study Site
Our experiment utilized a former agricultural field (~1 acre) at the University of Missouri’s Bradford Research Center in Columbia, MO (38.893604, -92.201154, Boone County, MO). Typical for the central U.S., the 10-year (2011 – 2021) mean annual precipitation and average air temperature at Bradford Research Center were 927.74 \(\pm\) 143.01 mm and 12.40 \(\pm\) 10.45 \(\circ\)C, respectively (Commercial Agriculture Automated Weather Station Network, d.n.). Prior to our experiment, the field we grew herbicide-resistant soybeans for at least three years, reflecting conditions similar to most prairie reconstructions before seeding (Rowe, 2010; Newbold et al., 2019). In March 2021, we tilled our study site and hand-removed rhizome clumps to create a smooth surface before our first seeding.
Experimental Design
We used seed additions to manipulate the arrival of 36 native tallgrass prairie plant species. Since prairies experience two peaks in dispersal activity (Rabinowitz & Rapp 1980), we classified species into two dispersal guilds, summer-dispersing species (first peak in dispersal activity before September 1st) and fall-dispersing species. Dispersal guild is more informative than flowering guild because dispersal does not always follow flowering (e.g., Penstemon digitalis). We based our classifications on expert opinion from the Missouri Department of Conservation, The Nature Conservancy, and our previous work on seed rain patterns in Missouri tallgrass prairies (Wynne et al. in prep). Species used in our study consisted of 29 native species captured in our previous study on seed rain and seven summer-dispersing species that restoration managers cited as having minimal success in prairie reconstructions (e.g., Viola pedatifida) (Table S1) (Barak et al., 2022). When possible, we obtained seeds from local ecotype commercial sellers. However, several summer-dispersing species were not grown commercially in Missouri and were sourced elsewhere. We stored the seeds in a refrigerator (2.78 \(\circ\)C) until seeding.
We used a randomized factorial block design to test the effects of seeding timing and order (Fig. 1). Each block (n = 6) contained five plots (2 x 2 m) randomly assigned a seed addition treatment following one of the four arrival treatments:
the addition of species in order of first peak in dispersal activity (NAT)
- May 24th 2021 (VIOPED, PACPLA) - June 4th 2021 (SISCAM, SPHOBT, CORLAN) - June 22nd 2021 (LOBSPI, TRAOHI, CARBUS, HEURIC) - July 11th 2021 (CORPAL, DODMEA) - July 25th 2021 (KOEMAC, AMOCAN) - August 2nd 2021 (LINSUL, MELVIR) - August 21st 2021 (DALCAN, DALPUR, ACHMIL) - September 17th 2021 (CROSAG) - October 4th 2021 (CHAFAS, MONFIS, RATPIN, SPOHET, RUDHIR) - October 18th 2021 (BIDARI, HELMOL, LESCAP) - November 1st 2021 (PENDIG, ERYYUC, HYPPUN, SORNUT, SCHSCO, LIAPYC) - November 19th 2021 (CORTRI, PYCTEN, SOLRIG)the lumped addition of 18 summer-dispersing species on March 22nd, 2021 followed by a lumped addition of 18 fall-dispersing species approximately five months later on September 5th, 2021 (LE)
the lumped addition of 18 fall-dispersing species on March 22nd, 2021 followed by a lumped addition of 18 summer-dispersing species approximately five months later on September 5th, 2021 (LL)
the simultaneous addition of the entire species pool on March 22nd, 2021 typical of most prairie reconstructions (SIM).
Additionally, we had an unseeded negative control (NON) to determine whether active seeding treatments improved community diversity compared to passive regeneration from ambient seed sources (e.g., seed rain and soil seed bank). We also seeded white clover (Trifolium repens) between the experimental plots to prevent erosion. In fall 2022, we started annually removing invading white clover and woody species from plots. All other unseeded species were not weeded to better reflect realistic assembly and reconstruction processes.
We started seeding at the end of the dormant season (March 22nd, 2021) and continued until November 19th, 2021. We hand-seeded species into experimental plots at a density of 50 seeds \(m^{-2}\) using sand as a broadcasting agent. For treatments requiring multiple seedings, we incorporated M-Binder tackifier (Ecology Controls, Carpinteria, CA), a natural adhesive often used in hydroseeding, into the seeding mixes to increase soil-seed contact without disturbing the existing vegetation. After seeding, we lightly watered plots to activate the tackifier. All plots received equal amounts of sand, tackifier, and water during every seeding.
Data Collection
Starting in 2022, we conducted floristic surveys to assess plant community diversity and composition in experimental plots for two years. To determine species abundance, we measured the percent aerial cover of all vascular plant species rooted within a 1 \(m^{2}\) subplot in the center of the experimental plot. Because many seeded summer-dispersing species senesce before peak biomass (August – September), we measured vegetative cover twice yearly, once in the early summer (late June) and again at peak biomass (late August). We only used the highest percent cover value for species present in both surveys. We identified species according to Yatskievych (1999, 2006, 2013). Additionally, we measured environmental variables including percent bare ground cover and light levels in plots. For each plot, we measured light levels using a (insert name of light bar) two times above the vegetation, in the midstory, and ground-level at solar noon \(\pm\) 2 hours.
Data Setup
Each tab below contains the code used to import, clean, and manipulate data to later use for further analyses.
Load libraries
## Libraries
### ----- Data Cleaning and Management -----
# Import excel files
library(readxl)
# Export dataframes to excel files
library(writexl)
# Data cleaning and visualization
library(tidyverse)
### ----- Data Analysis -----
# Community ecology functions (NMDS, diversity indices, etc.)
library(vegan)
library(BiodiversityR)
# Mixed models
library(lme4)
library(lmerTest)
library(MuMIn)
# Checks model assumptions
library(rstatix)
### ----- Data visualization -----
#Functions that assist in making high quality figures
library(cowplot)
library(ggpubr)
library(ggrepel)
library(flextable)Import Datasets
## Import Datasets
## inner (1 m^2) cover for Summer
#### comprehensive vegetative cover (%) for all species in the inner (1 m^2) permanent sampling area.
inner_cover_summer_2022 <- read_excel("~/Priority-Effects/Priority Effects - Github/Data/Cover/PE_Cover_Inner_Summer_2022_Cleaned.xlsx")
inner_cover_summer_2023 <- read_excel("~/Priority-Effects/Priority Effects - Github/Data/Cover/PE_Cover_Inner_Summer_2023_Cleaned.xlsx")
## inner (1 m^2) cover for Fall
#### comprehensive vegetative cover (%) for all species in the inner (1 m^2) permanent sampling area.
inner_cover_fall_2021 <- read_excel("~/Priority-Effects/Priority Effects - Github/Data/Cover/PE_Cover_Inner_Fall_2021_Cleaned.xlsx")
inner_cover_fall_2022 <- read_excel("~/Priority-Effects/Priority Effects - Github/Data/Cover/PE_Cover_Inner_Fall_2022_Cleaned.xlsx")
inner_cover_fall_2023 <- read_excel("~/Priority-Effects/Priority Effects - Github/Data/Cover/PE_Cover_Inner_Fall_2023_Cleaned.xlsx")
## Bradford weather
weather <- read_excel("~/Priority-Effects/Priority Effects - Github/Data/Environment/Bradford_10_Year_Weather.xlsx")
## Light
light_fall_2022 <- read_excel("~/Priority-Effects/Priority Effects - Github/Data/Environment/Priority_Effects_Fall_2022_Light.xlsx")
light_fall_2023 <- read_excel("~/Priority-Effects/Priority Effects - Github/Data/Environment/Priority_Effects_Fall_2023_Light.xlsx")
## Bare and litter - 2021
bare_cover_2021 <- read_excel("~/Priority-Effects/Priority Effects - Github/Data/Cover/PE_Bare_Ground_Cleaned_2021.xlsx")
## Species list
species_list_PE <- read_excel("~/Priority-Effects/Priority Effects - Github/Data/Species_List.xlsx")Dataset cleaning
2021
## 2021 Inner Cover - Dataset cleaning
### Remove unnecessary columns
#### Get rid of the notes, unknown, and senensced columns
### Fall
inner_cover_fall_2021 <- inner_cover_fall_2021[,-c(6,7,8)]
#### change log
# - senesced PLASPP changed to PLAVIR since PLALAN and PLAMAJ do not senesce until Oct. (Sep 6, 2023)
# - JUNSPP changed to JUNINT since JUNINT has been observed in 2022 and 2023 (Sep 6, 2023)
### Get rid of EMP
### Fall
inner_cover_fall_2021 <- subset(inner_cover_fall_2021, Treatment != "EMP")2022
## 2022 Inner Cover - Dataset cleaning
### Remove unnecessary columns
### Get rid of the notes, unknown, and sentenced columns
#### Summer
inner_cover_summer_2022 <- inner_cover_summer_2022[,-8]
#### change log
# - Fixed IDs: most GALAPA -> GALPED and HEURIC -> GEUCAN (Sep 6, 2023)
#### Fall
inner_cover_fall_2022 <- inner_cover_fall_2022[,-c(8,9,10)]
#### change log
# - Fixed IDs: VERSPP -> VERARV since this is the species observed in the summer (Sep 6, 2023)
### Get rid of EMP
#### Summer
inner_cover_summer_2022 <- subset(inner_cover_summer_2022, Treatment != "EMP")
#### Fall
inner_cover_fall_2022 <- subset(inner_cover_fall_2022, Treatment != "EMP")2023
## 2023 Inner Cover - Dataset cleaning
### Remove unnecessary columns
### Get rid of the notes, unknown, and sentenced columns
#### Summer
inner_cover_summer_2023 <- inner_cover_summer_2023[,-c(6,7)]
#### change log
# - Fixed IDs: was able to ID ALLSPP as ALLVIN (Sep 6, 2023)
#### Fall
inner_cover_fall_2023 <- inner_cover_fall_2023[,-c(6,7)]
### Get rid of EMP
#### Summer
inner_cover_summer_2023 <- subset(inner_cover_summer_2023, Treatment != "EMP")
#### Fall
inner_cover_fall_2023 <- subset(inner_cover_fall_2023, Treatment != "EMP")Joining Years
#Checking to make sure all the SPP6 levels are correct (no misspellings or weirdness)
### Fall 2021
sort(unique(inner_cover_fall_2021$SPP6))## [1] "ACAVIR" "ACHMIL" "AGRHYE" "AMATUB" "AMBART"
## [6] "AMBTRI" "BARVUL" "BIDARI" "BROJAP" "CARFRA"
## [11] "CHAFAS" "CONCAN" "CORLAN" "CORTIN" "CORTRI"
## [16] "CROSAG" "CYNDAC" "CYPODO" "DALCAN" "DIGISC"
## [21] "DIGSAN" "ECHSPP" "ERACIL" "ERIANN" "ERYYUC"
## [26] "GERCAR" "GLYMAX" "HELAUT" "HELMOL" "HORPUS"
## [31] "IPOHED" "IPOLAC" "IVAANN" "JUNINT" "LEPSPP"
## [36] "LESCAP" "LINSUL" "LONMAC" "LOTCOR" "MELOFF"
## [41] "MOLVER" "MONFIS" "MYOVER" "OXADIL" "PACKERA?"
## [46] "PANCAP" "PANDIC" "PENDIG" "PERHYD/PUN" "PERPEN"
## [51] "PHAARU" "PHYSPP" "PLALAN" "PLAMAJ" "PLAVIR"
## [56] "POTNOR" "RANARV" "RATPIN" "RUDHIR" "RUMCRI"
## [61] "SCHSCO" "SETFAB" "SETPUM" "SIDSPI" "SOLCAN"
## [66] "SOLCAR" "SORNUT" "SYMPIL" "TAROFF" "TRIREP"
## [71] "VERURT" "VIOPED" "XANSTR"### Summer 2022
sort(unique(inner_cover_summer_2022$SPP6))## [1] "ACHMIL" "AGRHYE" "BARVUL" "BROJAP" "CARBRE" "CARFRA" "CERSPP" "CORLAN"
## [9] "CORPAL" "ERIANN" "GALAPA" "GALPED" "GERCAR" "GEUCAN" "HORPUS" "LOTCOR"
## [17] "MELOFF" "MYOVER" "PENDIG" "PHAARU" "PLALAN" "PLAMAJ" "PLAVIR" "POTNOR"
## [25] "RANABO" "RUDHIR" "RUMCRI" "SILANT" "SPHOBT" "TAROFF" "TRIPRA" "TRIREP"
## [33] "VERARV" "VIOPED"### Fall 2022
sort(unique(inner_cover_fall_2022$SPP6))## [1] "ACAVIR" "ACHMIL" "AGRHYE" "AMBART" "AMBTRI"
## [6] "BARE" "BARVUL" "BIDARI" "BROJAP" "CARFRA"
## [11] "CARSPP" "CHAFAS" "CONCAN" "CORLAN" "CORPAL"
## [16] "CORTRI" "CYNDAC" "DESSPP" "DIGISC" "DIGSAN"
## [21] "ECHSPP" "ERASPE" "ERIANN" "ERYYUC" "FESARU"
## [26] "HELAUT" "HELMOL" "HORPUS" "IPOLAC" "IVAANN"
## [31] "JUNINT" "LACSER" "LESCAP" "LIAPYC" "LINSUL"
## [36] "LITTER" "LONMAC" "LOTCOR" "MONFIS" "MORALB"
## [41] "MYOVER" "OXADIL" "PANCAP" "PANDIC" "PENDIG"
## [46] "PERHYD/PUN" "PHAARU" "PLALAN" "PLAMAJ" "POAPRA"
## [51] "POTNOR" "RANABO" "RATPIN" "RUDHIR" "RUMCRI"
## [56] "SETFAB" "SETPUM" "SIDSPI" "SOLCAN" "SOLCAR"
## [61] "SORNUT" "SPOHET" "TAROFF" "TRIPRA" "TRIREP"
## [66] "VERARV" "VERURT" "VIOPED"### Summer 2023
sort(unique(inner_cover_summer_2023$SPP6))## [1] "ACHMIL" "AGRHYE" "ALLVIN" "BARVUL" "BROJAP" "CARBRE" "CARFRA" "CARSPP"
## [9] "CERSPP" "CORLAN" "ERIANN" "FESARU" "GALAPA" "GALPED" "GERCAR" "GEUCAN"
## [17] "HORPUS" "JUNINT" "KOEMAC" "LEPVIR" "LONMAC" "LOTCOR" "MELSPP" "MYOVER"
## [25] "OXADIL" "PENDIG" "PHAARU" "PLALAN" "PLAMAJ" "PLAVIR" "POAPRA" "PYCTEN"
## [33] "RANABO" "RUDHIR" "RUMCRI" "SILANT" "SPHOBT" "TAROFF" "TRIPRA" "TRIREP"
## [41] "VERARV" "VIOPED"### Fall 2023
sort(unique(inner_cover_fall_2023$SPP6))## [1] "ACHMIL" "AGRHYE" "AMBART" "AMBTRI" "BARE" "BARVUL" "BIDARI" "CARFRA"
## [9] "CARSPP" "CERSPP" "CHAFAS" "CONCAN" "CORLAN" "CORTRI" "CYNDAC" "DESSPP"
## [17] "ECHSPP" "ERASPE" "ERIANN" "ERYYUC" "FESARU" "GEUCAN" "HELAUT" "HELMOL"
## [25] "IPOHED" "IPOLAC" "IVAANN" "KOEMAC" "LESCAP" "LIAPYC" "LITTER" "LONMAC"
## [33] "LOTCOR" "MELALB" "MELSPP" "MONFIS" "MYOVER" "OXADIL" "PANDIC" "PENDIG"
## [41] "PHAARU" "PLALAN" "PLAMAJ" "POAPRA" "PYCTEN" "RANABO" "RATPIN" "RUDHIR"
## [49] "RUMCRI" "SETFAB" "SETPUM" "SIDSPI" "SOLCAN" "SOLCAR" "SORNUT" "SPOHET"
## [57] "SYMPIL" "TAROFF" "TRIPRA" "TRIREP" "VERURT" "VIOPED"### ^ above looks good so far# Make a unique identifier for year and season
### 2021
#### Season
inner_cover_fall_2021$Season <- rep("Fall", nrow=(inner_cover_fall_2021))
#### Year
inner_cover_fall_2021$Year <- rep("2021", nrow=(inner_cover_fall_2021))
### 2022
# Already has the identifiers
inner_cover_fall_2022$Year <- as.factor(inner_cover_fall_2022$Year)
inner_cover_summer_2022$Year <- as.factor(inner_cover_summer_2022$Year)
### 2023
#### Season
inner_cover_summer_2023$Season <- rep("Summer", nrow=(inner_cover_summer_2023))
inner_cover_fall_2023$Season <- rep("Fall", nrow=(inner_cover_fall_2023))
#### Year
inner_cover_fall_2023$Year <- rep("2023", nrow=(inner_cover_fall_2023))
inner_cover_fall_2023$Year <- as.factor(inner_cover_fall_2023$Year)
inner_cover_summer_2023$Year <- rep("2023", nrow=(inner_cover_summer_2023))
inner_cover_summer_2023$Year <- as.factor(inner_cover_summer_2023$Year)# Join summer and fall cover for each year
full_2021_data <- inner_cover_fall_2021
full_2022_data <- full_join(inner_cover_summer_2022,inner_cover_fall_2022)
full_2023_data <- full_join(inner_cover_summer_2023,inner_cover_fall_2023)
# Join all years together
### First 2021 and 2022
full_21_22_data <- full_join(full_2021_data, full_2022_data)
### Then 2021 and 2022 with 2023
full_year_data <- full_join(full_21_22_data, full_2023_data)# lump certain taxa together
### Lump the Carex for now (*****Ask Lauren about what to do about when the CARSPP > CARFRA or CARBRE?***)
## MELOFF + MELALB -> MELSPP
## LEPVIR -> LEPSPP
## CARFRA -> CARSPP
## CARBRE -> CARSPP
for(i in 1:nrow(full_year_data)) {
if(full_year_data[i,3] == "CARFRA"){full_year_data [i,3] <- "CARSPP"}
if(full_year_data[i,3] == "CARBRE"){full_year_data [i,3] <- "CARSPP"}
if(full_year_data[i,3] == "MELOFF"){full_year_data [i,3] <- "MELSPP"}
if(full_year_data[i,3] == "MELALB"){full_year_data [i,3] <- "MELSPP"}
if(full_year_data[i,3] == "LEPSPP"){full_year_data [i,3] <- "LEPVIR"}
if(full_year_data[i,3] == "ECHSPP"){full_year_data [i,3] <- "ECHCRU"}
# Make all cover data that was less than 1 = to 1 (to make the data discrete)
if(full_year_data[i,4] < 1) {full_year_data [i,4] <- 1}
}# Function below takes that highest cover value for a species seen twice (e.g., once in the summer and again in fall)
inner_cover_max <- full_year_data %>%
group_by(Block, Treatment, Year, SPP6) %>%
summarise(max_cover=max(Percent_Cover))
##### Note I manually checked species that were found in both the datasets (ACHMIL, TRIREP, CORLAN) and confirmed that the code above took the highest cover value for a species seen twice
### Remove Bare
inner_cover_max_only <- subset(inner_cover_max, SPP6 != "BARE")
### Remove Litter
inner_cover_max_only <- subset(inner_cover_max_only, SPP6 != "LITTER")#Checking again to make sure all the SPP6 levels are correct (no misspellings or weirdness)
sort(unique(inner_cover_max_only$SPP6))## [1] "ACAVIR" "ACHMIL" "AGRHYE" "ALLVIN" "AMATUB"
## [6] "AMBART" "AMBTRI" "BARVUL" "BIDARI" "BROJAP"
## [11] "CARSPP" "CERSPP" "CHAFAS" "CONCAN" "CORLAN"
## [16] "CORPAL" "CORTIN" "CORTRI" "CROSAG" "CYNDAC"
## [21] "CYPODO" "DALCAN" "DESSPP" "DIGISC" "DIGSAN"
## [26] "ECHCRU" "ERACIL" "ERASPE" "ERIANN" "ERYYUC"
## [31] "FESARU" "GALAPA" "GALPED" "GERCAR" "GEUCAN"
## [36] "GLYMAX" "HELAUT" "HELMOL" "HORPUS" "IPOHED"
## [41] "IPOLAC" "IVAANN" "JUNINT" "KOEMAC" "LACSER"
## [46] "LEPVIR" "LESCAP" "LIAPYC" "LINSUL" "LONMAC"
## [51] "LOTCOR" "MELSPP" "MOLVER" "MONFIS" "MORALB"
## [56] "MYOVER" "OXADIL" "PACKERA?" "PANCAP" "PANDIC"
## [61] "PENDIG" "PERHYD/PUN" "PERPEN" "PHAARU" "PHYSPP"
## [66] "PLALAN" "PLAMAJ" "PLAVIR" "POAPRA" "POTNOR"
## [71] "PYCTEN" "RANABO" "RANARV" "RATPIN" "RUDHIR"
## [76] "RUMCRI" "SCHSCO" "SETFAB" "SETPUM" "SIDSPI"
## [81] "SILANT" "SOLCAN" "SOLCAR" "SORNUT" "SPHOBT"
## [86] "SPOHET" "SYMPIL" "TAROFF" "TRIPRA" "TRIREP"
## [91] "VERARV" "VERURT" "VIOPED" "XANSTR"SPP6_list <- sort(unique(inner_cover_max_only$SPP6))
SPP6_list <- as.data.frame(SPP6_list)
### Saw at least 94 different species across study
#### More since I had to lump MELOFF, MELALB, CARFRA, and CARBRE
### Function below makes the species list
# write_xlsx(SPP6_list, "Species_List2.xlsx")### Add species information to the dataset
#### N/A is N = native, A = non-native, G = Genus
inner_cover_max_only <- full_join(inner_cover_max_only, species_list_PE)
#Remove unnecessary columns
inner_cover_max_reduced <- inner_cover_max_only[, -c(6,7,10)]Bare/Litter
### Bare dataset
inner_bare_cover <- subset(inner_cover_max, SPP6 == "BARE")
### Litter dataset
inner_litter_cover <- subset(inner_cover_max, SPP6 == "LITTER")
### Manipulate 2021 dataset to match 2022 and 2023
#### Bare Ground
bare_2021 <- subset(bare_cover_2021, Cover_Type != "Litter")
bare_2021 <- subset(bare_2021, Treatment != "EMP")
names(bare_2021) <- c("Block", "Treatment", "SPP6", "max_cover")
bare_2021$Year <- rep("2021", nrow(bare_2021))
for(i in 1:nrow(bare_2021)) {
if(bare_2021[i,3] == "Bare"){bare_2021[i,3] <- "BARE"}
}
#### Litter
litter_2021 <- subset(bare_cover_2021 , Cover_Type != "Bare")
litter_2021 <- subset(litter_2021, Treatment != "EMP")
names(litter_2021) <- c("Block", "Treatment", "SPP6", "max_cover")
litter_2021$Year <- rep("2021", nrow(litter_2021))
for(i in 1:nrow(litter_2021)) {
if(bare_2021[i,3] == "Litter"){bare_2021[i,3] <- "LITTER"}
}
### Full_join datasets
inner_bare_cover <- full_join(inner_bare_cover, bare_2021)
inner_litter_cover <- full_join(inner_litter_cover, litter_2021)Question 1: Does timing of species arrival influence local diversity?
Setup and Data Management
Data Analysis
We fit mixed-effects linear models with block as a random effect to determine whether year, timing of species arrival via seeding treatments, and the interaction between year and seeding treatments influenced local diversity (i.e., species richness, mean conservatism value, and Shannon diversity index) in reconstructed tallgrass prairie plant communities. In cases where random effect variance was estimated as near zero, we dropped the random effect and fit a simpler linear model instead. For each diversity model, we conducted a type III analysis of variance (ANOVA) followed by a posthoc Tukey test to determine significant pairwise differences.
Below is the code used to calculate the species richness, mean conservatism value, and Shannon diversity index for each plot, treatment, and year combination.
### Making a summary table of diversity indices between treatments
#Make a unique identifier for every block and treatment
inner_cover_max_reduced$UniqueBlock <- paste(inner_cover_max_reduced$Treatment, inner_cover_max_reduced$Block, sep="_")
### Making a community matrix
inner_cover_max_reduced_lumped <- inner_cover_max_reduced %>%
filter(Year != "2021") %>%
group_by(Block, Treatment, Year, SPP6) %>%
summarise(max_cover = sum(max_cover))
#Make a wide formatted dataset
inner_wide <- inner_cover_max_reduced_lumped %>%
spread(key="SPP6", value="max_cover")
#Replace NA values with 0
inner_wide[is.na(inner_wide)] <- 0
#Make a separate dataframe for the labels
inner_wide.labs <- inner_wide[, c(1,2,3)]
#Turn our dataset into a matrix
inner_wide_mat <- inner_wide[,-c(1,2,3)]
#Double check your work using the View() function.
#View(inner_wide_mat)## Calculate diversity metrics
# --- Species richness ---
Species_richness <- specnumber(inner_wide_mat)
# --- Shannon's Diversity ---
#Calculate Shannon's Diversity
Shannon <- diversity(inner_wide_mat, index="shannon")
# --- Mean C ---
# Make a dataset that has each treatment, the species found in that treatment, and their associated C values
inner_cover_max_reduced$C_Value <- as.numeric(inner_cover_max_reduced$C_Value)
C_hist <- inner_cover_max_reduced %>%
filter( C_Value >= 0) %>%
group_by(Block, Treatment, SPP6, Year, C_Value) %>%
summarize(tot.cover = sum(max_cover)) %>%
select(-c(tot.cover))
# Calculate mean C for each Treatment in a Block
### Remember you have to filter out anything that isn't native!
Only_C <- inner_cover_max_reduced %>%
filter(Year != "2021") %>%
filter( C_Value >= 0) %>%
filter( C_Value != "NA")
Only_C$C_Value <- as.numeric(Only_C$C_Value)
Mean_C <- Only_C %>%
group_by(Block, Treatment, Year) %>%
summarize(Mean_C = mean(C_Value))Diversity Indices - Summary Statistics
Below is a summary table showing the mean and SD for mean C value (C), species richness (SR), and Shannon diversity index (SDI) for each seeding treatment.
Treatment | Year | C | ± SD | Richness | ± SD | SDI | ± SD |
|---|---|---|---|---|---|---|---|
LE | 2022 | 1.68 | 0.58 | 20.00 | 3.46 | 2.12 | 0.20 |
LE | 2023 | 1.91 | 0.47 | 18.17 | 1.83 | 2.18 | 0.16 |
LL | 2022 | 1.99 | 0.41 | 21.00 | 2.19 | 2.25 | 0.20 |
LL | 2023 | 2.62 | 0.20 | 17.83 | 2.40 | 2.16 | 0.13 |
NAT | 2022 | 1.43 | 0.84 | 17.67 | 2.34 | 2.05 | 0.12 |
NAT | 2023 | 1.45 | 0.73 | 17.00 | 3.69 | 2.01 | 0.28 |
NON | 2022 | 0.48 | 0.34 | 15.00 | 2.61 | 1.86 | 0.17 |
NON | 2023 | 0.42 | 0.41 | 14.50 | 2.74 | 1.75 | 0.20 |
SIM | 2022 | 2.34 | 0.44 | 24.50 | 3.89 | 2.23 | 0.32 |
SIM | 2023 | 2.84 | 0.32 | 22.17 | 5.08 | 2.34 | 0.36 |
Mean C
BioR.theme <- theme(
panel.background = element_blank(),
panel.border = element_blank(),
panel.grid = element_blank(),
axis.line = element_line("gray25"),
text = element_text(size = 16, family="Arial"),
axis.text = element_text(size = 13, colour = "gray25"),
axis.title = element_text(size = 16, colour = "gray25"),
legend.title = element_text(size = 16),
legend.text = element_text(size = 16),
legend.key = element_blank())Mean C - Model assumptions
Distribution of C values for all treatments were right skewed, indicating the majority of native plants in all treatments were of low conservatism value. Simultaneous seeding had the most “even” distribution of C values and almost every value was represented in this treatment. Across years, the largest changes in the distribution of C values were in the treatments that seeded the fall dispersing species in September. These changes were attributed to the replacement of species with low c values with more conservative species. Plots that were never seeded mainly consisted of species with low conservatism value.
Generally, mean C values were normally distributed with few extreme outliers, suggesting we can use a linear model to predict mean c as a function of treatment and year.
Mean C - Model fitting
I used a linear model to predict mean c as a function of seeding treatment, year, and an interaction between treatment and year. A mixed model using block as a random effect produced a singular fit because the variance explained by the random effect block was estimated as close to zero and did not further inform the data. Therefore, I chose to use the simpler model.
Model fit was high, with an adjusted \(R^{2}\) value of 0.687. Also, residual plots did not reveal any concerning patterns.
##
## Call:
## lm(formula = Mean_C ~ Treatment + Year, data = Summary_table)
##
## Residuals:
## Min 1Q Median 3Q Max
## -1.07124 -0.32383 -0.07948 0.28266 1.12107
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 1.6602 0.1608 10.327 2.16e-14 ***
## TreatmentLL 0.5149 0.2075 2.481 0.016251 *
## TreatmentNAT -0.3549 0.2075 -1.710 0.093003 .
## TreatmentNON -1.3445 0.2075 -6.479 2.89e-08 ***
## TreatmentSIM 0.7983 0.2075 3.847 0.000318 ***
## Year2023 0.2660 0.1313 2.026 0.047691 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.5084 on 54 degrees of freedom
## Multiple R-squared: 0.7134, Adjusted R-squared: 0.6869
## F-statistic: 26.88 on 5 and 54 DF, p-value: 1.544e-13# Residual plot looks decent
plot(Mean_C.mod.lm.additive)
Analysis of variance (type III) revealed that seeding treatment was the only significant factor influence mean c value.
## Anova Table (Type III tests)
##
## Response: Mean_C
## Sum Sq Df F value Pr(>F)
## (Intercept) 27.561 1 106.6537 2.160e-14 ***
## Treatment 33.674 4 32.5772 8.088e-14 ***
## Year 1.061 1 4.1057 0.04769 *
## Residuals 13.955 54
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Posthoc analysis indicated that all seeding treatments increased mean c value compared to unseeded plots. Seeding all species in one addition produced more conservative communities than treatments that seeded summer dispersing species first followed by fall dispersing species. Even though seeding only the fall dispersing species first produced higher quality communities than those that seeded species with natural timing, there was no difference in mean c between the seeding treatments that had two lumped additions.
Mean_C.mod.lm.est <- lm(Mean_C~Treatment, data =Summary_table)
est.lm.treatment_treatment <- emmeans::emmeans(Mean_C.mod.lm.est, ~ Treatment, type = "response")
pairs(est.lm.treatment_treatment, adjust = "tukey")## contrast estimate SE df t.ratio p.value
## LE - LL -0.515 0.213 55 -2.414 0.1270
## LE - NAT 0.355 0.213 55 1.664 0.4644
## LE - NON 1.345 0.213 55 6.303 <.0001
## LE - SIM -0.798 0.213 55 -3.742 0.0039
## LL - NAT 0.870 0.213 55 4.077 0.0013
## LL - NON 1.859 0.213 55 8.717 <.0001
## LL - SIM -0.283 0.213 55 -1.329 0.6748
## NAT - NON 0.990 0.213 55 4.639 0.0002
## NAT - SIM -1.153 0.213 55 -5.406 <.0001
## NON - SIM -2.143 0.213 55 -10.045 <.0001
##
## P value adjustment: tukey method for comparing a family of 5 estimates#Posthoc test
#est.lm.treatment <- emmeans::emmeans(Mean_C.mod.lm.interaction, ~ Treatment*Year, type = "response")
#pairs(est.lm.treatment, adjust = "tukey")Mean C - Figures
Overall, seeding fall dispersing species early either as part of one simultaneous or as two lumped additions produced the highest quality communities.
Species Richness
SIM treatments had the greatest species richness compared to all other treatments, while NON and NAT had the lowest.
SR - Model assumptions
Overall, species richness follows a normal distribution despite being count data. There were a couple of outliers in the natural seeding treatment but this is likely alright.
SR - Model fitting
I chose to use a linear mixed-effects model with block as a random effect to predict species richness as a function of seeding treatment, year, and an interaction between treatment and year. Despite being count data, the distribution of species richness is remarkably balanced and not skewed. Model comparison using AIC revealed that the Normal model performed better than the Poisson model, which was underdispersed. A more complicated Quasipoisson model that accounted for underdispersion produced similar results to the Normal mixed model. Therefore, I chose to use the less complicated Normal model.
Overall, model fit was decent. For example, \(R^{2}_{GLMM(m)}\) (i.e., the variance explained by only the fixed effects) was 0.474 and \(R^{2}_{GLMM(c)}\) (i.e., variance explained by the entire model) was 0.598. Inclusion of the random effect block substantially improved model fit. Furthermore, there were no strange patterns or occurrences in the model residuals.
## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: Species_richness ~ Treatment + Year + (1 | Block)
## Data: Summary_table
##
## REML criterion at convergence: 231.5
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -3.0095 -0.7513 0.1272 0.6940 1.5317
##
## Random effects:
## Groups Name Variance Std.Dev.
## Block (Intercept) 1.828 1.352
## Residual 8.513 2.918
## Number of obs: 48, groups: Block, 6
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) 18.3333 1.0915 23.7757 16.796 1.09e-14 ***
## TreatmentLE 1.7500 1.1912 38.0000 1.469 0.1500
## TreatmentLL 2.0833 1.1912 38.0000 1.749 0.0884 .
## TreatmentSIM 6.0000 1.1912 38.0000 5.037 1.18e-05 ***
## Year2023 -2.0000 0.8423 38.0000 -2.375 0.0227 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) TrtmLE TrtmLL TrtSIM
## TreatmentLE -0.546
## TreatmentLL -0.546 0.500
## TreatmntSIM -0.546 0.500 0.500
## Year2023 -0.386 0.000 0.000 0.000## R2m R2c
## [1,] 0.3644925 0.4768066
A type III ANOVA revealed significant differences in species richness between years and seeding treatments. However, there was no significant interactive effect between year and seeding treatment on species richness.
Seeding treatments that added species either in one or two lumped additions had significantly greater species richness than unseeded plots. Multiple seedings according to natural dispersal activity did not increase species richness compared to unseeded plots. Adding the entire species pool simultaneously produced communities with the greatest species richness. Treatments that conducted multiple seedings resulted in communities possessing similar amounts of species richness regardless of timing and order of species arrival.
## Analysis of Deviance Table (Type III Wald chisquare tests)
##
## Response: Species_richness
## Chisq Df Pr(>Chisq)
## (Intercept) 282.1205 1 < 2.2e-16 ***
## Treatment 27.1051 3 5.596e-06 ***
## Year 5.6383 1 0.01757 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1## contrast estimate SE df t.ratio p.value
## NAT - LE -1.750 1.26 39 -1.389 0.5137
## NAT - LL -2.083 1.26 39 -1.653 0.3616
## NAT - SIM -6.000 1.26 39 -4.762 0.0002
## LE - LL -0.333 1.26 39 -0.265 0.9934
## LE - SIM -4.250 1.26 39 -3.373 0.0088
## LL - SIM -3.917 1.26 39 -3.108 0.0177
##
## Degrees-of-freedom method: kenward-roger
## P value adjustment: tukey method for comparing a family of 4 estimatesSR - Figures
Overall, seeding all species at the same time produced communities with the greatest species richness. For the most part, seeding increased species richness
Summer-Dispersing
### Filter out everything that isn't a seeded summer-dispersing species
inner_long <- inner_wide %>%
gather(key = "SPP6", value = "max_cover", ACAVIR:VIOPED)
Sown_species_ID <- left_join(inner_long, species_list_PE)## Joining with `by = join_by(SPP6)`SR_summer <- Sown_species_ID %>%
group_by(Block, Treatment, Year) %>%
filter(Seeded == "Summer" & max_cover > 0 & Treatment != "NON") %>%
summarize(Summer_SR = n_distinct(SPP6))## `summarise()` has grouped output by 'Block', 'Treatment'. You can override
## using the `.groups` argument.SR_summer_zero <- Sown_species_ID %>%
group_by(Block, Treatment, Year) %>%
filter(Seeded == "Summer" & Treatment != "NON") %>%
summarise(max_cover = sum(max_cover))## `summarise()` has grouped output by 'Block', 'Treatment'. You can override
## using the `.groups` argument.SR_summer_zero <- SR_summer_zero %>%
filter(max_cover == 0)
SR_summer_zero <- SR_summer_zero [,c(1,2,3)]
SR_summer_zero$Summer_SR <- rep(0, nrow(SR_summer_zero))
SR_summer_all <- full_join(SR_summer,SR_summer_zero)## Joining with `by = join_by(Block, Treatment, Year, Summer_SR)`## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: Summer_SR ~ Treatment + Year + (1 | Block)
## Data: SR_summer_all
##
## REML criterion at convergence: 133.4
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -2.3183 -0.5433 -0.2825 0.5178 2.1606
##
## Random effects:
## Groups Name Variance Std.Dev.
## Block (Intercept) 0.06321 0.2514
## Residual 0.92763 0.9631
## Number of obs: 48, groups: Block, 6
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) 2.2917 0.3274 34.2203 7.000 4.31e-08 ***
## TreatmentLL -1.9167 0.3932 38.0000 -4.875 1.96e-05 ***
## TreatmentNAT -1.9167 0.3932 38.0000 -4.875 1.96e-05 ***
## TreatmentSIM -0.1667 0.3932 38.0000 -0.424 0.674
## Year2023 0.2500 0.2780 38.0000 0.899 0.374
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) TrtmLL TrtNAT TrtSIM
## TreatmentLL -0.601
## TreatmntNAT -0.601 0.500
## TreatmntSIM -0.601 0.500 0.500
## Year2023 -0.425 0.000 0.000 0.000## R2m R2c
## [1,] 0.4697125 0.5035431Anova(SR.summer.mod.normal, type = 3)## Analysis of Deviance Table (Type III Wald chisquare tests)
##
## Response: Summer_SR
## Chisq Df Pr(>Chisq)
## (Intercept) 49.0067 1 2.551e-12 ***
## Treatment 43.6596 3 1.783e-09 ***
## Year 0.8085 1 0.3686
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1SR.summer.mod.normal.est <- lmer(Summer_SR~Treatment+(1|Block), data =
SR_summer_all)
est.lmer.SR.summer.treatment <- emmeans::emmeans(SR.summer.mod.normal.est, ~ Treatment, type = "response")
pairs(est.lmer.SR.summer.treatment , adjust = "tukey")## contrast estimate SE df t.ratio p.value
## LE - LL 1.917 0.392 39 4.887 0.0001
## LE - NAT 1.917 0.392 39 4.887 0.0001
## LE - SIM 0.167 0.392 39 0.425 0.9739
## LL - NAT 0.000 0.392 39 0.000 1.0000
## LL - SIM -1.750 0.392 39 -4.462 0.0004
## NAT - SIM -1.750 0.392 39 -4.462 0.0004
##
## Degrees-of-freedom method: kenward-roger
## P value adjustment: tukey method for comparing a family of 4 estimatesset.seed(7)
SR_summer_all_mean_sd <- SR_summer_all %>%
group_by(Treatment, Year) %>%
summarize(mean_SR_summer = mean(Summer_SR), sd_SR_summer = sd(Summer_SR))## `summarise()` has grouped output by 'Treatment'. You can override using the
## `.groups` argument.SR_summer_all_mean_sd$Treatment <- factor(SR_summer_all_mean_sd$Treatment, levels=c("NAT", "LE", "LL", "SIM"))
SR_summer_all$Treatment <- factor(SR_summer_all$Treatment, levels=c("NAT", "LE", "LL", "SIM"))
sr.summer.plot_treatment <- ggplot()+
geom_jitter(data = SR_summer_all,
aes(x = Year, y = Summer_SR, fill = Treatment,
shape = Treatment, group = Treatment),
alpha = 1, size = 3.5,
position = position_jitterdodge(dodge.width = 0.85))+
geom_errorbar(data = SR_summer_all_mean_sd,
aes(x = Year, y = mean_SR_summer, group = Treatment,
ymin = mean_SR_summer - sd_SR_summer,
ymax = mean_SR_summer + sd_SR_summer),
position = position_dodge(width = 0.85), width = .5,
size = 0.75)+
geom_point(data = SR_summer_all_mean_sd, size = 4.5,
aes(x = Year, y = mean_SR_summer,
shape = Treatment, group = Treatment),
position = position_dodge2(width = 0.85),
fill = "black")+
labs(y = "Sown summer richness", x = "" )+
theme_classic() +
theme(text=element_text(size=18),
legend.key.size=unit(1, "cm"),
legend.position="none",
axis.text.x = element_text(size = 16),
axis.text.y = element_text(size = 16))+
scale_fill_manual(name = "Treatment",
breaks = c("LE", "LL", "NAT", "SIM"),
labels = c("Summer first", "Fall first", "Natural", "Simultaneous"),
values = c("#009E73", "#56B4E9", "#F0E442", "#0072B2"))+
# scale_x_discrete(breaks = c("LE", "LL", "NAT", "SIM"),
# labels = c("Summer first", "Fall first", "Natural", "Simultaneous"))+
scale_shape_manual(name = "Treatment",
breaks = c("LE", "LL", "NAT", "SIM"),
labels = c("Summer first", "Fall first", "Natural", "Simultaneous"),
values=c(21, 22,23, 25))+
scale_y_continuous(breaks = seq(0, 20, by = 5),
limits = c(-0.75, 20))
sr.summer.plot_treatment <- sr.summer.plot_treatment +
annotate("text", x = 1.2, y = 20, label = "Treatment: p < 0.001***",
fontface = 2, size = 4)Fall-Dispersing
### Filter out everything that isn't a seeded fall-dispersing species
SR_fall <- Sown_species_ID %>%
filter(Seeded == "Fall" & max_cover > 0 & Treatment != "NON") %>%
summarize(Fall_SR = n_distinct(SPP6))## `summarise()` has grouped output by 'Block', 'Treatment'. You can override
## using the `.groups` argument.SR_fall_zero <- Sown_species_ID %>%
group_by(Block, Treatment, Year) %>%
filter(Seeded == "Fall" & Treatment != "NON") %>%
summarise(max_cover = sum(max_cover))## `summarise()` has grouped output by 'Block', 'Treatment'. You can override
## using the `.groups` argument.SR_fall_zero <- SR_fall_zero %>%
filter(max_cover == 0)
SR_fall_zero <- SR_fall_zero[,c(1,2,3)]
SR_fall_zero$Fall_SR <- rep(0, nrow(SR_fall_zero))
SR_fall_all <- full_join(SR_fall,SR_fall_zero)## Joining with `by = join_by(Block, Treatment, Year, Fall_SR)`##
## Call:
## lm(formula = Fall_SR ~ Treatment + Year, data = SR_fall_all)
##
## Residuals:
## Min 1Q Median 3Q Max
## -2.9583 -0.7083 0.0417 0.7292 3.3750
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 2.7083 0.4977 5.442 2.36e-06 ***
## TreatmentLL 3.8333 0.6295 6.089 2.71e-07 ***
## TreatmentNAT 0.5000 0.6295 0.794 0.431
## TreatmentSIM 5.1667 0.6295 8.207 2.42e-10 ***
## Year2023 -0.2500 0.4452 -0.562 0.577
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 1.542 on 43 degrees of freedom
## Multiple R-squared: 0.6919, Adjusted R-squared: 0.6633
## F-statistic: 24.15 on 4 and 43 DF, p-value: 1.608e-10Anova(SR.fall.mod.normal , type = 3)## Anova Table (Type III tests)
##
## Response: Fall_SR
## Sum Sq Df F value Pr(>F)
## (Intercept) 70.417 1 29.6129 2.356e-06 ***
## Treatment 228.917 3 32.0894 4.753e-11 ***
## Year 0.750 1 0.3154 0.5773
## Residuals 102.250 43
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1SR.fall.mod.normal.est <- lm(Fall_SR~Treatment, data = SR_fall_all)
est.lmer.SR.fall.treatment <- emmeans::emmeans(SR.fall.mod.normal.est, ~ Treatment, type = "response")
pairs(est.lmer.SR.fall.treatment , adjust = "tukey")## contrast estimate SE df t.ratio p.value
## LE - LL -3.83 0.625 44 -6.137 <.0001
## LE - NAT -0.50 0.625 44 -0.800 0.8538
## LE - SIM -5.17 0.625 44 -8.272 <.0001
## LL - NAT 3.33 0.625 44 5.337 <.0001
## LL - SIM -1.33 0.625 44 -2.135 0.1582
## NAT - SIM -4.67 0.625 44 -7.471 <.0001
##
## P value adjustment: tukey method for comparing a family of 4 estimatesset.seed(5)
SR_fall_all_mean_sd <- SR_fall_all %>%
group_by(Treatment, Year) %>%
summarize(mean_SR_fall = mean(Fall_SR), sd_SR_fall = sd(Fall_SR))## `summarise()` has grouped output by 'Treatment'. You can override using the
## `.groups` argument.SR_fall_all_mean_sd $Treatment <- factor(SR_fall_all_mean_sd $Treatment, levels=c("NAT", "LE", "LL", "SIM"))
SR_fall_all$Treatment <- factor(SR_fall_all$Treatment, levels=c("NAT", "LE", "LL", "SIM"))
sr.fall.plot_treatment <- ggplot()+
geom_jitter(data = SR_fall_all,
aes(x = Year, y = Fall_SR, fill = Treatment, shape = Treatment, group = Treatment),
alpha = 0.8, size = 3.5,
position = position_jitterdodge(dodge.width = 0.85))+
geom_errorbar(data = SR_fall_all_mean_sd,
aes(x = Year, y = mean_SR_fall, ymin = mean_SR_fall - sd_SR_fall,
ymax = mean_SR_fall + sd_SR_fall, group = Treatment),
position = position_dodge(width = 0.85), width = .5, size = 0.75) +
geom_point(data = SR_fall_all_mean_sd, size = 4.5,
aes(x = Year, y = mean_SR_fall, shape = Treatment, group = Treatment),
fill = "black", position = position_dodge(width = 0.85))+
labs(y = "Sown fall richness", x = "" ) +
theme_classic() +
theme(text=element_text(size=18),
legend.key.size=unit(1, "cm"),
legend.position="none",
axis.text.x = element_text(size = 16),
axis.text.y = element_text(size = 16))+
scale_fill_manual(name = "Treatment",
breaks = c("LE", "LL", "NAT", "SIM"),
labels = c("Summer first", "Fall first", "Natural", "Simultaneous"),
values = c("#009E73", "#56B4E9", "#F0E442", "#0072B2"))+
scale_shape_manual(name = "Treatment",
breaks = c("LE", "LL", "NAT", "SIM"),
labels = c("Summer first", "Fall first", "Natural", "Simultaneous"),
values=c(21, 22,23, 25))+
scale_y_continuous(breaks = seq(0, 20, by = 5),
limits = c(-0.75, 20))
sr.fall.plot_treatment <- sr.fall.plot_treatment +
annotate("text", x = 1.2, y = 20, label = "Treatment: p < 0.001***",
fontface = 2, size = 4)Unsown
### Filter out everything that isn't a seeded fall-dispersing species
SR_unsown_all <- Sown_species_ID %>%
filter(Seeded == "Unseeded" & max_cover > 0 & Treatment != "NON") %>%
summarize(unsown_SR = n_distinct(SPP6))## `summarise()` has grouped output by 'Block', 'Treatment'. You can override
## using the `.groups` argument.## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: unsown_SR ~ Treatment + Year + (1 | Block)
## Data: SR_unsown_all
##
## REML criterion at convergence: 212.9
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -2.11194 -0.60308 0.06422 0.64409 1.79134
##
## Random effects:
## Groups Name Variance Std.Dev.
## Block (Intercept) 0.5141 0.717
## Residual 5.8202 2.413
## Number of obs: 48, groups: Block, 6
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) 15.0833 0.8318 32.4490 18.132 < 2e-16 ***
## TreatmentLL -1.5833 0.9849 38.0000 -1.608 0.11620
## TreatmentNAT -0.3333 0.9849 38.0000 -0.338 0.73689
## TreatmentSIM -0.7500 0.9849 38.0000 -0.761 0.45106
## Year2023 -2.0000 0.6964 38.0000 -2.872 0.00664 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) TrtmLL TrtNAT TrtSIM
## TreatmentLL -0.592
## TreatmntNAT -0.592 0.500
## TreatmntSIM -0.592 0.500 0.500
## Year2023 -0.419 0.000 0.000 0.000Anova(SR.unsown.mod.normal , type = 3)## Analysis of Deviance Table (Type III Wald chisquare tests)
##
## Response: unsown_SR
## Chisq Df Pr(>Chisq)
## (Intercept) 328.7867 1 < 2.2e-16 ***
## Treatment 2.8922 3 0.408540
## Year 8.2472 1 0.004082 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1set.seed(5)
SR_unsown_all_mean_sd <- SR_unsown_all %>%
group_by(Treatment, Year) %>%
summarize(mean_SR_unsown = mean(unsown_SR), sd_SR_unsown = sd(unsown_SR))## `summarise()` has grouped output by 'Treatment'. You can override using the
## `.groups` argument.SR_unsown_all_mean_sd $Treatment <- factor(SR_unsown_all_mean_sd $Treatment, levels=c("NAT", "LE", "LL", "SIM"))
SR_unsown_all$Treatment <- factor(SR_unsown_all$Treatment, levels=c("NAT", "LE", "LL", "SIM"))
sr.unsown.plot_treatment <- ggplot()+
geom_jitter(data = SR_unsown_all,
aes(x = Year, y = unsown_SR, fill = Treatment, shape = Treatment, group = Treatment),
alpha = 0.8, size = 3.5,
position = position_jitterdodge(dodge.width = 0.85))+
geom_errorbar(data = SR_unsown_all_mean_sd,
aes(x = Year, y = mean_SR_unsown, ymin = mean_SR_unsown - sd_SR_unsown,
ymax = mean_SR_unsown + sd_SR_unsown, group = Treatment),
position = position_dodge(width = 0.85), width = .5, size = 0.75) +
geom_point(data = SR_unsown_all_mean_sd, size = 4.5,
aes(x = Year, y = mean_SR_unsown, shape = Treatment, group = Treatment),
fill = "black", position = position_dodge(width = 0.85))+
labs(y = "Unsown richness", x = "" ) +
theme_classic() +
theme(text=element_text(size=18),
legend.key.size=unit(1, "cm"),
legend.position="none",
axis.text.x = element_text(size = 16),
axis.text.y = element_text(size = 16))+
scale_fill_manual(name = "Treatment",
breaks = c("LE", "LL", "NAT", "SIM"),
labels = c("Summer first", "Fall first", "Natural", "Simultaneous"),
values = c("#009E73", "#56B4E9", "#F0E442", "#0072B2"))+
scale_shape_manual(name = "Treatment",
breaks = c("LE", "LL", "NAT", "SIM"),
labels = c("Summer first", "Fall first", "Natural", "Simultaneous"),
values=c(21, 22,23, 25))+
scale_y_continuous(breaks = seq(0, 20, by = 5),
limits = c(-0.75, 20))
sr.unsown.plot_treatment
Shannon Diversity Index
SDI - Model assumptions
SIM and LL treatments had the greatest Shannon diversity index value while NON had the lowest. NAT and LE appear comparable.
SDI - Model fitting
## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: Shannon ~ Treatment + Year + (1 | Block)
## Data: Summary_table
##
## REML criterion at convergence: 7.8
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -2.42601 -0.49701 0.08675 0.60144 1.94713
##
## Random effects:
## Groups Name Variance Std.Dev.
## Block (Intercept) 0.004618 0.06796
## Residual 0.049227 0.22187
## Number of obs: 48, groups: Block, 6
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) 2.026960 0.076796 31.983342 26.394 < 2e-16 ***
## TreatmentLE 0.122187 0.090579 38.000000 1.349 0.18533
## TreatmentLL 0.173774 0.090579 38.000000 1.918 0.06259 .
## TreatmentSIM 0.251361 0.090579 38.000000 2.775 0.00851 **
## Year2023 0.008802 0.064049 38.000000 0.137 0.89142
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) TrtmLE TrtmLL TrtSIM
## TreatmentLE -0.590
## TreatmentLL -0.590 0.500
## TreatmntSIM -0.590 0.500 0.500
## Year2023 -0.417 0.000 0.000 0.000## Analysis of Deviance Table (Type III Wald chisquare tests)
##
## Response: Shannon
## Chisq Df Pr(>Chisq)
## (Intercept) 696.6555 1 < 2e-16 ***
## Treatment 8.1465 3 0.04308 *
## Year 0.0189 1 0.89070
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1## R2m R2c
## [1,] 0.1370613 0.2110774## contrast estimate SE df t.ratio p.value
## NAT - LE -0.1222 0.0906 38 -1.349 0.5384
## NAT - LL -0.1738 0.0906 38 -1.918 0.2375
## NAT - SIM -0.2514 0.0906 38 -2.775 0.0406
## LE - LL -0.0516 0.0906 38 -0.570 0.9405
## LE - SIM -0.1292 0.0906 38 -1.426 0.4913
## LL - SIM -0.0776 0.0906 38 -0.857 0.8269
##
## Results are averaged over the levels of: Year
## Degrees-of-freedom method: kenward-roger
## P value adjustment: tukey method for comparing a family of 4 estimates
SDI - Figures
Plot Panel
Question 2: Does timing of species arrival influence seeded cover?
note some of the species used in the seed mix were present in the seed bank at my study site (e.g., RUDHIR and PENDIG). Additionally, some spillover occurred with seeded species seeding into other treatments like NON (e.g., BIDARI)
Seeded2 <- Seeded %>% filter(Year == "2023") %>% filter(max_cover > 0)
unique(Seeded2$SPP6)## [1] "ACHMIL" "BIDARI" "CHAFAS" "CORLAN" "CORTRI" "ERYYUC" "HELMOL" "KOEMAC"
## [9] "LESCAP" "LIAPYC" "MONFIS" "PENDIG" "PYCTEN" "RATPIN" "RUDHIR" "SORNUT"
## [17] "SPHOBT" "SPOHET" "VIOPED"Total Seeded Cover
TSC - Model fitting
##
## Call:
## glm(formula = cbind(seeded_cover, tot_cover) ~ Treatment * Year,
## family = binomial, data = Total_Cover_Tot)
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -2.11998 0.08881 -23.870 < 2e-16 ***
## TreatmentLL 0.94324 0.10412 9.059 < 2e-16 ***
## TreatmentNAT -0.38112 0.13799 -2.762 0.00575 **
## TreatmentSIM 1.06270 0.10397 10.221 < 2e-16 ***
## Year2023 0.62783 0.11282 5.565 2.63e-08 ***
## TreatmentLL:Year2023 -0.41302 0.13685 -3.018 0.00254 **
## TreatmentNAT:Year2023 -0.36533 0.18174 -2.010 0.04441 *
## TreatmentSIM:Year2023 -0.37777 0.13662 -2.765 0.00569 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for binomial family taken to be 1)
##
## Null deviance: 837.79 on 47 degrees of freedom
## Residual deviance: 293.23 on 40 degrees of freedom
## AIC: 557.58
##
## Number of Fisher Scoring iterations: 5## Analysis of Deviance Table (Type III tests)
##
## Response: cbind(seeded_cover, tot_cover)
## LR Chisq Df Pr(>Chisq)
## Treatment 272.445 3 < 2.2e-16 ***
## Year 32.035 1 1.514e-08 ***
## Treatment:Year 10.305 3 0.01614 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1## contrast odds.ratio SE df null z.ratio p.value
## LE Year2022 / LL Year2022 0.389 0.0405 Inf 1 -9.059 <.0001
## LE Year2022 / NAT Year2022 1.464 0.2020 Inf 1 2.762 0.1048
## LE Year2022 / SIM Year2022 0.346 0.0359 Inf 1 -10.221 <.0001
## LE Year2022 / LE Year2023 0.534 0.0602 Inf 1 -5.565 <.0001
## LE Year2022 / LL Year2023 0.314 0.0328 Inf 1 -11.076 <.0001
## LE Year2022 / NAT Year2023 1.126 0.1469 Inf 1 0.909 0.9853
## LE Year2022 / SIM Year2023 0.269 0.0281 Inf 1 -12.574 <.0001
## LL Year2022 / NAT Year2022 3.760 0.4466 Inf 1 11.150 <.0001
## LL Year2022 / SIM Year2022 0.887 0.0680 Inf 1 -1.558 0.7751
## LL Year2022 / LE Year2023 1.371 0.1210 Inf 1 3.572 0.0085
## LL Year2022 / LL Year2023 0.807 0.0625 Inf 1 -2.774 0.1016
## LL Year2022 / NAT Year2023 2.892 0.3181 Inf 1 9.654 <.0001
## LL Year2022 / SIM Year2023 0.691 0.0534 Inf 1 -4.784 <.0001
## NAT Year2022 / SIM Year2022 0.236 0.0280 Inf 1 -12.169 <.0001
## NAT Year2022 / LE Year2023 0.365 0.0461 Inf 1 -7.977 <.0001
## NAT Year2022 / LL Year2023 0.215 0.0256 Inf 1 -12.917 <.0001
## NAT Year2022 / NAT Year2023 0.769 0.1096 Inf 1 -1.842 0.5909
## NAT Year2022 / SIM Year2023 0.184 0.0219 Inf 1 -14.231 <.0001
## SIM Year2022 / LE Year2023 1.545 0.1361 Inf 1 4.935 <.0001
## SIM Year2022 / LL Year2023 0.909 0.0702 Inf 1 -1.234 0.9218
## SIM Year2022 / NAT Year2023 3.259 0.3580 Inf 1 10.753 <.0001
## SIM Year2022 / SIM Year2023 0.779 0.0600 Inf 1 -3.246 0.0258
## LE Year2023 / LL Year2023 0.588 0.0523 Inf 1 -5.971 <.0001
## LE Year2023 / NAT Year2023 2.109 0.2495 Inf 1 6.312 <.0001
## LE Year2023 / SIM Year2023 0.504 0.0447 Inf 1 -7.729 <.0001
## LL Year2023 / NAT Year2023 3.585 0.3958 Inf 1 11.563 <.0001
## LL Year2023 / SIM Year2023 0.857 0.0667 Inf 1 -1.988 0.4901
## NAT Year2023 / SIM Year2023 0.239 0.0264 Inf 1 -12.981 <.0001
##
## P value adjustment: tukey method for comparing a family of 8 estimates
## Tests are performed on the log odds ratio scale## contrast odds.ratio SE df null z.ratio p.value
## LE / LL 0.479 0.0328 Inf 1 -10.767 <.0001
## LE / NAT 1.757 0.1597 Inf 1 6.204 <.0001
## LE / SIM 0.417 0.0285 Inf 1 -12.792 <.0001
## LL / NAT 3.671 0.2977 Inf 1 16.039 <.0001
## LL / SIM 0.872 0.0476 Inf 1 -2.510 0.0584
## NAT / SIM 0.237 0.0192 Inf 1 -17.751 <.0001
##
## Results are averaged over the levels of: Year
## P value adjustment: tukey method for comparing a family of 4 estimates
## Tests are performed on the log odds ratio scale## contrast odds.ratio SE df null z.ratio p.value
## Year2022 / Year2023 0.713 0.0378 Inf 1 -6.391 <.0001
##
## Results are averaged over the levels of: Treatment
## Tests are performed on the log odds ratio scale
TSC - Figures
Seeded cover of summer dispersing species
ESC - Model fitting
##
## Call:
## glm(formula = cbind(seeded_cover, tot_cover) ~ Treatment + Year,
## family = binomial, data = Total_Cover_Early)
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -2.9526 0.1056 -27.955 < 2e-16 ***
## TreatmentLL -3.3561 0.3411 -9.838 < 2e-16 ***
## TreatmentNAT -2.0984 0.1995 -10.519 < 2e-16 ***
## TreatmentSIM -0.3623 0.1070 -3.387 0.000707 ***
## Year2023 1.0821 0.1106 9.780 < 2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for binomial family taken to be 1)
##
## Null deviance: 734.86 on 47 degrees of freedom
## Residual deviance: 247.79 on 43 degrees of freedom
## AIC: 381.05
##
## Number of Fisher Scoring iterations: 6## Analysis of Deviance Table (Type III tests)
##
## Response: cbind(seeded_cover, tot_cover)
## LR Chisq Df Pr(>Chisq)
## Treatment 373.50 3 < 2.2e-16 ***
## Year 106.69 1 < 2.2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1## contrast odds.ratio SE df null z.ratio p.value
## LE / LL 28.6759 9.7820 Inf 1 9.838 <.0001
## LE / NAT 8.1535 1.6265 Inf 1 10.519 <.0001
## LE / SIM 1.4366 0.1537 Inf 1 3.387 0.0039
## LL / NAT 0.2843 0.1089 Inf 1 -3.284 0.0056
## LL / SIM 0.0501 0.0172 Inf 1 -8.706 <.0001
## NAT / SIM 0.1762 0.0360 Inf 1 -8.502 <.0001
##
## Results are averaged over the levels of: Year
## P value adjustment: tukey method for comparing a family of 4 estimates
## Tests are performed on the log odds ratio scale## contrast odds.ratio SE df null z.ratio p.value
## Year2022 / Year2023 0.339 0.0375 Inf 1 -9.780 <.0001
##
## Results are averaged over the levels of: Treatment
## Tests are performed on the log odds ratio scale
ESC - Figures
Fall dispersing species only
LSC - Model fitting
## Generalized linear mixed model fit by maximum likelihood (Laplace
## Approximation) [glmerMod]
## Family: binomial ( logit )
## Formula: cbind(seeded_cover, tot_cover) ~ Treatment + Year + (1 | Block)
## Data: Total_Cover_Late
##
## AIC BIC logLik deviance df.resid
## 605.3 616.5 -296.6 593.3 42
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -4.3147 -1.7714 -0.6057 1.2574 6.5425
##
## Random effects:
## Groups Name Variance Std.Dev.
## Block (Intercept) 0.008129 0.09016
## Number of obs: 48, groups: Block, 6
##
## Fixed effects:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -2.74863 0.09364 -29.352 < 2e-16 ***
## TreatmentLL 1.59358 0.09078 17.554 < 2e-16 ***
## TreatmentNAT 0.16919 0.11161 1.516 0.12953
## TreatmentSIM 1.55200 0.09191 16.885 < 2e-16 ***
## Year2023 0.13908 0.05056 2.751 0.00594 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) TrtmLL TrtNAT TrtSIM
## TreatmentLL -0.796
## TreatmntNAT -0.643 0.664
## TreatmntSIM -0.789 0.807 0.656
## Year2023 -0.279 0.018 -0.002 0.024## Analysis of Deviance Table (Type III Wald chisquare tests)
##
## Response: cbind(seeded_cover, tot_cover)
## Chisq Df Pr(>Chisq)
## (Intercept) 861.5598 1 < 2e-16 ***
## Treatment 566.5093 3 < 2e-16 ***
## Year 7.5683 1 0.00594 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1## contrast odds.ratio SE df null z.ratio p.value
## LE / LL 0.203 0.0184 Inf 1 -17.554 <.0001
## LE / NAT 0.844 0.0942 Inf 1 -1.516 0.4278
## LE / SIM 0.212 0.0195 Inf 1 -16.885 <.0001
## LL / NAT 4.155 0.3536 Inf 1 16.738 <.0001
## LL / SIM 1.042 0.0592 Inf 1 0.732 0.8841
## NAT / SIM 0.251 0.0216 Inf 1 -16.031 <.0001
##
## Results are averaged over the levels of: Year
## P value adjustment: tukey method for comparing a family of 4 estimates
## Tests are performed on the log odds ratio scale## contrast odds.ratio SE df null z.ratio p.value
## Year2022 / Year2023 0.87 0.044 Inf 1 -2.751 0.0059
##
## Results are averaged over the levels of: Treatment
## Tests are performed on the log odds ratio scale
LSC - Figures
Plot Panel
Species specific differences
##
## One Sample t-test
##
## data: test2$difference
## t = 5.4921, df = 89, p-value = 3.722e-07
## alternative hypothesis: true mean is not equal to 0
## 95 percent confidence interval:
## 3.935637 8.397697
## sample estimates:
## mean of x
## 6.166667
Question 3: Does timing of species arrival influence species composition? - Bray-Curtis
Year: 2022
## Call: capscale(formula = inner.hel.2022 ~ Treatment, data =
## env.data.2022, distance = "bray", add = TRUE)
##
## Inertia Proportion Rank
## Total 3.7536 1.0000
## Constrained 0.8657 0.2306 3
## Unconstrained 2.8879 0.7694 20
## Inertia is Lingoes adjusted squared Bray distance
## Species scores projected from 'inner.hel.2022'
##
## Eigenvalues for constrained axes:
## CAP1 CAP2 CAP3
## 0.5388 0.2425 0.0843
##
## Eigenvalues for unconstrained axes:
## MDS1 MDS2 MDS3 MDS4 MDS5 MDS6 MDS7 MDS8
## 0.5403 0.4052 0.2960 0.2529 0.1945 0.1649 0.1568 0.1338
## (Showing 8 of 20 unconstrained eigenvalues)
##
## Constant added to distances: 0.03513658## Permutation test for capscale under reduced model
## Terms added sequentially (first to last)
## Permutation: free
## Number of permutations: 999
##
## Model: capscale(formula = inner.hel.2022 ~ Treatment, data = env.data.2022, distance = "bray", add = TRUE)
## Df SumOfSqs F Pr(>F)
## Treatment 3 0.86567 1.9984 0.001 ***
## Residual 20 2.88789
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1## Multiple comparisons for capscale for all contrasts of Treatment
##
## Model: BiodiversityR::multiconstrained(method = "capscale", formula = inner.hel.2022 ~ Treatment, data = env.data.2022, distance = "bray", add = TRUE, by = "term")
##
## Df SumOfSqs F Pr(>F)
## LE vs. LL 1 0.32432 3.2161 0.002 **
## LE vs. NAT 1 0.08466 0.6738 0.812
## LE vs. SIM 1 0.31377 2.6708 0.007 **
## LL vs. NAT 1 0.25386 2.4115 0.006 **
## LL vs. SIM 1 0.17717 1.9078 0.030 *
## NAT vs. SIM 1 0.37810 3.0297 0.006 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## UniqueBlock Block Treatment Year Bare_cover Above Ground Middle axis1
## 1 LE_1 1 LE <NA> NA 1575.5 67 239.5 -0.07771696
## 2 LL_1 1 LL 2022 14 1565.0 31 215.0 0.48336017
## 3 NAT_1 1 NAT 2022 4 1535.0 32 128.5 -0.43746584
## 4 SIM_1 1 SIM 2022 2 1580.5 22 626.0 0.58299824
## 5 LE_2 2 LE 2022 10 1611.0 163 422.0 -0.61581200
## 6 LL_2 2 LL 2022 11 1601.5 71 336.0 0.27390902
## axis2 labels
## 1 0.38755734 LE_1
## 2 -0.95497148 LL_1
## 3 -0.03453055 NAT_1
## 4 1.97668133 SIM_1
## 5 1.09725065 LE_2
## 6 -0.60130836 LL_2## axis1 axis2 labels
## ACAVIR -0.013759760 0.0105436448 ACAVIR
## ACHMIL -0.039905108 0.3151005315 ACHMIL
## AGRHYE 0.008944552 0.0420245595 AGRHYE
## AMBART -0.165234521 0.0893240283 AMBART
## AMBTRI -0.395292407 -0.4196119068 AMBTRI
## BARVUL -0.022532498 -0.0940021827 BARVUL
## BIDARI 0.616701687 -0.0381609416 BIDARI
## BROJAP 0.113385176 0.0873793317 BROJAP
## CARSPP -0.133319955 -0.0709646462 CARSPP
## CERSPP 0.013460740 0.0070992909 CERSPP
## CHAFAS 0.221973715 -0.0615936378 CHAFAS
## CONCAN 0.026314578 0.0186473255 CONCAN
## CORLAN 0.047431342 0.1387802756 CORLAN
## CORPAL -0.043650430 0.0230115005 CORPAL
## CORTRI 0.259290736 0.0956570002 CORTRI
## DESSPP -0.002361261 -0.0276468238 DESSPP
## DIGISC -0.005805795 0.0523668346 DIGISC
## DIGSAN 0.004596870 0.0292496220 DIGSAN
## ECHCRU -0.007543745 -0.0139571869 ECHCRU
## ERASPE 0.017765474 0.0161270608 ERASPE
## ERIANN -0.020170112 -0.0512415405 ERIANN
## ERYYUC 0.034816004 -0.0109188717 ERYYUC
## FESARU 0.001119958 -0.0145731064 FESARU
## GALPED -0.008890104 0.0004354018 GALPED
## GERCAR -0.029033113 -0.0210649548 GERCAR
## GEUCAN -0.018517585 -0.0096680542 GEUCAN
## HELAUT -0.105661417 -0.0149522445 HELAUT
## HELMOL 0.129128316 -0.0683750610 HELMOL
## HORPUS 0.035115845 -0.0733476653 HORPUS
## IPOLAC -0.006344513 0.0278887241 IPOLAC
## IVAANN 0.262803154 -0.0578486576 IVAANN
## JUNINT 0.011148351 -0.0176837652 JUNINT
## LACSER 0.017301947 0.0157062823 LACSER
## LESCAP 0.116802974 -0.0116205743 LESCAP
## LIAPYC 0.017301947 0.0157062823 LIAPYC
## LINSUL -0.013759760 0.0105436448 LINSUL
## LONMAC 0.004421824 0.0270484397 LONMAC
## LOTCOR -0.082750521 0.3856035839 LOTCOR
## MELSPP -0.043445692 -0.0226830500 MELSPP
## MONFIS -0.009543532 0.0202951310 MONFIS
## MYOVER -0.039536260 0.0611374149 MYOVER
## OXADIL -0.032682355 0.0243796838 OXADIL
## PANCAP 0.003201855 0.0279832615 PANCAP
## PANDIC -0.003845438 -0.0200239744 PANDIC
## PENDIG -0.082854822 0.0077086447 PENDIG
## PERHYD/PUN -0.015929940 0.0122065818 PERHYD/PUN
## PHAARU -0.106205632 -0.0554500935 PHAARU
## PLALAN 0.037094780 -0.0381481848 PLALAN
## PLAMAJ -0.037542985 0.0089781276 PLAMAJ
## PLAVIR -0.036194605 -0.0025266870 PLAVIR
## POAPRA 0.046096825 0.0669222606 POAPRA
## POTNOR 0.034618018 -0.0117608580 POTNOR
## RANABO 0.063310648 0.0623983669 RANABO
## RATPIN 0.295359023 -0.1903300280 RATPIN
## RUDHIR 0.128747186 -0.1095167813 RUDHIR
## RUMCRI -0.128181852 0.0250011764 RUMCRI
## SETFAB -0.087935392 -0.0525684638 SETFAB
## SETPUM -0.027218251 0.0588510004 SETPUM
## SIDSPI 0.019576599 0.0177711555 SIDSPI
## SOLCAN -0.015360372 -0.0080196692 SOLCAN
## SOLCAR 0.013348165 -0.0229292736 SOLCAR
## SORNUT 0.423351229 -0.0195015817 SORNUT
## SPHOBT 0.031491428 0.0285871434 SPHOBT
## SPOHET 0.025712643 0.0233413048 SPOHET
## TAROFF 0.005036684 0.0524681788 TAROFF
## TRIPRA 0.040655258 0.0369058434 TRIPRA
## TRIREP -0.379489398 -0.0327028056 TRIREP
## VERARV -0.035216468 0.0164490440 VERARV
## VERURT 0.019309511 -0.0306291798 VERURT
## VIOPED -0.017935072 0.0434831986 VIOPED## axis ggplot label
## 1 1 xlab.label CAP1 (14.4%)
## 2 2 ylab.label CAP2 (6.5%)## r p axis1 axis2 labels
## ACAVIR 0.01156079 0.916 -0.013759760 0.0105436448 ACAVIR
## ACHMIL 0.26416995 0.036 -0.039905108 0.3151005315 ACHMIL
## AGRHYE 0.04221726 0.640 0.008944552 0.0420245595 AGRHYE
## AMBART 0.05133255 0.574 -0.165234521 0.0893240283 AMBART
## AMBTRI 0.62167784 0.001 -0.395292407 -0.4196119068 AMBTRI
## BARVUL 0.24118220 0.062 -0.022532498 -0.0940021827 BARVUL
## BIDARI 0.80347542 0.001 0.616701687 -0.0381609416 BIDARI
## BROJAP 0.11728601 0.306 0.113385176 0.0873793317 BROJAP
## CARSPP 0.08699372 0.381 -0.133319955 -0.0709646462 CARSPP
## CERSPP 0.01879381 0.846 0.013460740 0.0070992909 CERSPP
## CHAFAS 0.32314850 0.016 0.221973715 -0.0615936378 CHAFAS
## CONCAN 0.05347984 0.563 0.026314578 0.0186473255 CONCAN
## CORLAN 0.23679300 0.068 0.047431342 0.1387802756 CORLAN
## CORPAL 0.22152803 0.072 -0.043650430 0.0230115005 CORPAL
## CORTRI 0.57092202 0.002 0.259290736 0.0956570002 CORTRI
## DESSPP 0.08205552 0.426 -0.002361261 -0.0276468238 DESSPP
## DIGISC 0.06806011 0.495 -0.005805795 0.0523668346 DIGISC
## DIGSAN 0.31708125 0.016 0.004596870 0.0292496220 DIGSAN
## ECHCRU 0.01043125 0.912 -0.007543745 -0.0139571869 ECHCRU
## ERASPE 0.16914332 0.090 0.017765474 0.0161270608 ERASPE
## ERIANN 0.10723478 0.337 -0.020170112 -0.0512415405 ERIANN
## ERYYUC 0.17864825 0.129 0.034816004 -0.0109188717 ERYYUC
## FESARU 0.01976619 0.825 0.001119958 -0.0145731064 FESARU
## GALPED 0.07257478 0.453 -0.008890104 0.0004354018 GALPED
## GERCAR 0.03530864 0.695 -0.029033113 -0.0210649548 GERCAR
## GEUCAN 0.07878386 0.578 -0.018517585 -0.0096680542 GEUCAN
## HELAUT 0.02234030 0.824 -0.105661417 -0.0149522445 HELAUT
## HELMOL 0.32842446 0.010 0.129128316 -0.0683750610 HELMOL
## HORPUS 0.09919497 0.349 0.035115845 -0.0733476653 HORPUS
## IPOLAC 0.09105548 0.401 -0.006344513 0.0278887241 IPOLAC
## IVAANN 0.44187949 0.004 0.262803154 -0.0578486576 IVAANN
## JUNINT 0.08948316 0.472 0.011148351 -0.0176837652 JUNINT
## LACSER 0.01396845 0.859 0.017301947 0.0157062823 LACSER
## LESCAP 0.24964718 0.045 0.116802974 -0.0116205743 LESCAP
## LIAPYC 0.01396845 0.859 0.017301947 0.0157062823 LIAPYC
## LINSUL 0.01156079 0.916 -0.013759760 0.0105436448 LINSUL
## LONMAC 0.09228958 0.362 0.004421824 0.0270484397 LONMAC
## LOTCOR 0.60180742 0.001 -0.082750521 0.3856035839 LOTCOR
## MELSPP 0.12993356 0.128 -0.043445692 -0.0226830500 MELSPP
## MONFIS 0.09716618 0.320 -0.009543532 0.0202951310 MONFIS
## MYOVER 0.10341326 0.323 -0.039536260 0.0611374149 MYOVER
## OXADIL 0.04738271 0.598 -0.032682355 0.0243796838 OXADIL
## PANCAP 0.13538690 0.231 0.003201855 0.0279832615 PANCAP
## PANDIC 0.04478159 0.623 -0.003845438 -0.0200239744 PANDIC
## PENDIG 0.09639392 0.349 -0.082854822 0.0077086447 PENDIG
## PERHYD/PUN 0.05224616 0.742 -0.015929940 0.0122065818 PERHYD/PUN
## PHAARU 0.01682139 0.855 -0.106205632 -0.0554500935 PHAARU
## PLALAN 0.07651040 0.445 0.037094780 -0.0381481848 PLALAN
## PLAMAJ 0.13775527 0.204 -0.037542985 0.0089781276 PLAMAJ
## PLAVIR 0.03838884 0.702 -0.036194605 -0.0025266870 PLAVIR
## POAPRA 0.36792330 0.006 0.046096825 0.0669222606 POAPRA
## POTNOR 0.02307529 0.791 0.034618018 -0.0117608580 POTNOR
## RANABO 0.36519063 0.007 0.063310648 0.0623983669 RANABO
## RATPIN 0.75052506 0.001 0.295359023 -0.1903300280 RATPIN
## RUDHIR 0.22610230 0.069 0.128747186 -0.1095167813 RUDHIR
## RUMCRI 0.14928772 0.185 -0.128181852 0.0250011764 RUMCRI
## SETFAB 0.06697542 0.477 -0.087935392 -0.0525684638 SETFAB
## SETPUM 0.08685742 0.384 -0.027218251 0.0588510004 SETPUM
## SIDSPI 0.30323356 0.048 0.019576599 0.0177711555 SIDSPI
## SOLCAN 0.12993356 0.128 -0.015360372 -0.0080196692 SOLCAN
## SOLCAR 0.02744986 0.774 0.013348165 -0.0229292736 SOLCAR
## SORNUT 0.70787270 0.001 0.423351229 -0.0195015817 SORNUT
## SPHOBT 0.11670926 0.333 0.031491428 0.0285871434 SPHOBT
## SPOHET 0.11670926 0.333 0.025712643 0.0233413048 SPOHET
## TAROFF 0.08199427 0.406 0.005036684 0.0524681788 TAROFF
## TRIPRA 0.11670926 0.333 0.040655258 0.0369058434 TRIPRA
## TRIREP 0.45984733 0.002 -0.379489398 -0.0327028056 TRIREP
## VERARV 0.03991363 0.650 -0.035216468 0.0164490440 VERARV
## VERURT 0.08948316 0.472 0.019309511 -0.0306291798 VERURT
## VIOPED 0.01282063 0.878 -0.017935072 0.0434831986 VIOPED## r p axis1 axis2 labels
## AMBTRI 0.6216778 0.001 -0.39529241 -0.41961191 AMBTRI
## BIDARI 0.8034754 0.001 0.61670169 -0.03816094 BIDARI
## CORTRI 0.5709220 0.002 0.25929074 0.09565700 CORTRI
## IVAANN 0.4418795 0.004 0.26280315 -0.05784866 IVAANN
## LOTCOR 0.6018074 0.001 -0.08275052 0.38560358 LOTCOR
## POAPRA 0.3679233 0.006 0.04609683 0.06692226 POAPRA
## RANABO 0.3651906 0.007 0.06331065 0.06239837 RANABO
## RATPIN 0.7505251 0.001 0.29535902 -0.19033003 RATPIN
## SORNUT 0.7078727 0.001 0.42335123 -0.01950158 SORNUT
## TRIREP 0.4598473 0.002 -0.37948940 -0.03270281 TRIREP
2022 NMDS
adonis2(inner.hel.2022 ~ Treatment, data = env.data.2022, dist="bray")## Permutation test for adonis under reduced model
## Terms added sequentially (first to last)
## Permutation: free
## Number of permutations: 999
##
## adonis2(formula = inner.hel.2022 ~ Treatment, data = env.data.2022, dist = "bray")
## Df SumOfSqs R2 F Pr(>F)
## Treatment 3 0.76026 0.25812 2.3195 0.001 ***
## Residual 20 2.18516 0.74188
## Total 23 2.94542 1.00000
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1### Pair-wise comparisons
library(pairwiseAdonis)## Loading required package: clusterpairwise.adonis2(inner.hel.2022 ~ Treatment, data = env.data.2022, p.adjust = "bonferroni")## $parent_call
## [1] "inner.hel.2022 ~ Treatment , strata = Null , permutations 999"
##
## $LE_vs_LL
## Df SumOfSqs R2 F Pr(>F)
## Treatment 1 0.32134 0.24719 3.2835 0.004 **
## Residual 10 0.97863 0.75281
## Total 11 1.29997 1.00000
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## $LE_vs_NAT
## Df SumOfSqs R2 F Pr(>F)
## Treatment 1 0.08466 0.06313 0.6738 0.79
## Residual 10 1.25652 0.93687
## Total 11 1.34118 1.00000
##
## $LE_vs_SIM
## Df SumOfSqs R2 F Pr(>F)
## Treatment 1 0.31377 0.21079 2.6708 0.012 *
## Residual 10 1.17481 0.78921
## Total 11 1.48859 1.00000
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## $LL_vs_NAT
## Df SumOfSqs R2 F Pr(>F)
## Treatment 1 0.24962 0.19812 2.4707 0.006 **
## Residual 10 1.01035 0.80188
## Total 11 1.25997 1.00000
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## $LL_vs_SIM
## Df SumOfSqs R2 F Pr(>F)
## Treatment 1 0.17717 0.16022 1.9078 0.037 *
## Residual 10 0.92864 0.83978
## Total 11 1.10581 1.00000
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## $NAT_vs_SIM
## Df SumOfSqs R2 F Pr(>F)
## Treatment 1 0.37395 0.23661 3.0994 0.009 **
## Residual 10 1.20653 0.76339
## Total 11 1.58048 1.00000
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## attr(,"class")
## [1] "pwadstrata" "list"Year: 2023
## Call: capscale(formula = inner.hel.2023 ~ Treatment, data =
## env.data.2023, distance = "bray", add = TRUE)
##
## Inertia Proportion Rank
## Total 4.9591 1.0000
## Constrained 1.3768 0.2776 3
## Unconstrained 3.5823 0.7224 20
## Inertia is Lingoes adjusted squared Bray distance
## Species scores projected from 'inner.hel.2023'
##
## Eigenvalues for constrained axes:
## CAP1 CAP2 CAP3
## 0.7696 0.4783 0.1289
##
## Eigenvalues for unconstrained axes:
## MDS1 MDS2 MDS3 MDS4 MDS5 MDS6 MDS7 MDS8
## 0.7364 0.4885 0.3927 0.3212 0.2513 0.2101 0.1807 0.1551
## (Showing 8 of 20 unconstrained eigenvalues)
##
## Constant added to distances: 0.05150478## Permutation test for capscale under reduced model
## Terms added sequentially (first to last)
## Permutation: free
## Number of permutations: 999
##
## Model: capscale(formula = inner.hel.2023 ~ Treatment, data = env.data.2023, distance = "bray", add = TRUE)
## Df SumOfSqs F Pr(>F)
## Treatment 3 1.3768 2.5622 0.001 ***
## Residual 20 3.5823
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1## Multiple comparisons for capscale for all contrasts of Treatment
##
## Model: BiodiversityR::multiconstrained(method = "capscale", formula = inner.hel.2023 ~ Treatment, data = env.data.2023, distance = "bray", add = TRUE)
##
## Df SumOfSqs F Pr(>F)
## LE vs. LL 1 0.56594 4.8862 0.002 **
## LE vs. NAT 1 0.22834 1.5823 0.098 .
## LE vs. SIM 1 0.45445 3.3335 0.003 **
## LL vs. NAT 1 0.34988 2.7638 0.004 **
## LL vs. SIM 1 0.28030 2.5274 0.004 **
## NAT vs. SIM 1 0.59997 3.6140 0.003 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## UniqueBlock Block Treatment Year Bare_cover Above Ground Middle axis1
## 1 LE_1 1 LE 2023 3 1852.0 114.0 1305.5 -0.1492502
## 2 LL_1 1 LL 2023 3 1840.0 59.0 1363.0 0.5871939
## 3 NAT_1 1 NAT 2023 2 1854.5 127.0 1335.0 -0.6337099
## 4 SIM_1 1 SIM 2023 13 1819.5 87.5 1486.5 1.0284389
## 5 LE_2 2 LE 2023 18 1786.0 196.5 1252.0 -0.8450437
## 6 LL_2 2 LL 2023 46 1731.0 109.5 248.0 0.3190052
## axis2 labels
## 1 0.8082363 LE_1
## 2 -1.1020888 LL_1
## 3 -0.6006288 NAT_1
## 4 0.9438200 SIM_1
## 5 0.8878608 LE_2
## 6 -0.4246796 LL_2## axis1 axis2 labels
## ACHMIL -0.0711105232 0.592520826 ACHMIL
## AGRHYE 0.0085664677 0.050758260 AGRHYE
## ALLVIN -0.0185623066 -0.016238688 ALLVIN
## AMBART -0.1280267347 0.032816374 AMBART
## AMBTRI -0.1574108999 -0.300615993 AMBTRI
## BARVUL -0.0269443322 -0.077481104 BARVUL
## BIDARI 0.0033968330 0.028002451 BIDARI
## BROJAP -0.0633261442 -0.103598384 BROJAP
## CARSPP -0.1212440666 -0.161028841 CARSPP
## CERSPP 0.0181368471 0.018000893 CERSPP
## CHAFAS 0.1606510797 0.055077603 CHAFAS
## CONCAN -0.0277846622 -0.052836330 CONCAN
## CORLAN -0.0212794604 0.241260701 CORLAN
## CORTRI 0.4820760468 0.112414607 CORTRI
## DESSPP -0.0157311613 -0.013651272 DESSPP
## ECHCRU 0.0223030327 0.019201888 ECHCRU
## ERASPE 0.0216923651 -0.019686356 ERASPE
## ERIANN -0.0132088157 -0.043483443 ERIANN
## ERYYUC 0.1372907701 0.002771190 ERYYUC
## FESARU -0.0135022460 -0.002617965 FESARU
## GALAPA -0.0146881661 0.013442529 GALAPA
## GALPED 0.0224825505 0.019236293 GALPED
## GERCAR 0.0108559078 -0.010491515 GERCAR
## GEUCAN -0.0142451724 -0.012361753 GEUCAN
## HELAUT -0.0583436417 -0.086514015 HELAUT
## HELMOL 0.3553255256 -0.034761103 HELMOL
## HORPUS -0.1475521122 -0.293480395 HORPUS
## IPOLAC -0.0090387556 -0.031331641 IPOLAC
## IVAANN 0.0455330518 0.039201868 IVAANN
## JUNINT 0.0136331062 -0.016777789 JUNINT
## KOEMAC -0.0487436827 0.170858654 KOEMAC
## LEPVIR -0.0572340552 0.003521346 LEPVIR
## LESCAP 0.1504633256 -0.021946339 LESCAP
## LIAPYC 0.0179812798 0.015481057 LIAPYC
## LONMAC 0.0004196733 0.031553343 LONMAC
## LOTCOR -0.0022111948 0.443190652 LOTCOR
## MONFIS -0.0331070975 0.064020687 MONFIS
## MYOVER -0.0378162987 0.046423163 MYOVER
## OXADIL -0.0721611624 0.020946220 OXADIL
## PANDIC -0.0207573701 0.008207711 PANDIC
## PENDIG -0.1192823538 0.025241468 PENDIG
## PHAARU -0.0859544746 -0.074776638 PHAARU
## PLALAN 0.0202676327 0.017207299 PLALAN
## PLAMAJ -0.0162510996 0.014872917 PLAMAJ
## PLAVIR -0.0141967833 0.012992818 PLAVIR
## POAPRA 0.0492069968 0.092727229 POAPRA
## PYCTEN -0.0189514955 -0.016445831 PYCTEN
## RANABO 0.0487820743 0.041907584 RANABO
## RATPIN 0.5240609418 -0.268029972 RATPIN
## RUDHIR -0.0093949298 -0.096825923 RUDHIR
## RUMCRI -0.0807045564 0.043057146 RUMCRI
## SETFAB -0.0327750423 -0.097865804 SETFAB
## SETPUM -0.4458985710 0.039650195 SETPUM
## SIDSPI -0.0161814285 -0.014042007 SIDSPI
## SOLCAR -0.0193228840 -0.027701254 SOLCAR
## SORNUT 0.7621760713 0.011031048 SORNUT
## SPHOBT 0.0372025018 0.032029647 SPHOBT
## SPOHET 0.0440449605 0.037920690 SPOHET
## SYMPIL -0.0346702003 0.027790858 SYMPIL
## TAROFF -0.0852401230 -0.004968236 TAROFF
## TRIPRA 0.0672797885 0.057924811 TRIPRA
## TRIREP -0.1045058986 0.087732099 TRIREP
## VERARV 0.0221452014 0.016543448 VERARV
## VERURT 0.0136331062 -0.016777789 VERURT
## VIOPED -0.0363283830 0.033247537 VIOPED## axis ggplot label
## 1 1 xlab.label CAP1 (15.5%)
## 2 2 ylab.label CAP2 (9.6%)## r p axis1 axis2 labels
## ACHMIL 0.557034573 0.001 -0.0711105232 0.592520826 ACHMIL
## AGRHYE 0.043162268 0.634 0.0085664677 0.050758260 AGRHYE
## ALLVIN 0.042659805 0.627 -0.0185623066 -0.016238688 ALLVIN
## AMBART 0.065102575 0.492 -0.1280267347 0.032816374 AMBART
## AMBTRI 0.332015883 0.018 -0.1574108999 -0.300615993 AMBTRI
## BARVUL 0.283979345 0.016 -0.0269443322 -0.077481104 BARVUL
## BIDARI 0.003289345 0.964 0.0033968330 0.028002451 BIDARI
## BROJAP 0.027646890 0.743 -0.0633261442 -0.103598384 BROJAP
## CARSPP 0.121420182 0.245 -0.1212440666 -0.161028841 CARSPP
## CERSPP 0.009249716 0.908 0.0181368471 0.018000893 CERSPP
## CHAFAS 0.287760954 0.019 0.1606510797 0.055077603 CHAFAS
## CONCAN 0.147106291 0.173 -0.0277846622 -0.052836330 CONCAN
## CORLAN 0.364442071 0.014 -0.0212794604 0.241260701 CORLAN
## CORTRI 0.635936293 0.001 0.4820760468 0.112414607 CORTRI
## DESSPP 0.027217602 0.918 -0.0157311613 -0.013651272 DESSPP
## ECHCRU 0.032565663 0.892 0.0223030327 0.019201888 ECHCRU
## ERASPE 0.045055339 0.611 0.0216923651 -0.019686356 ERASPE
## ERIANN 0.084905162 0.402 -0.0132088157 -0.043483443 ERIANN
## ERYYUC 0.211941307 0.092 0.1372907701 0.002771190 ERYYUC
## FESARU 0.017652522 0.812 -0.0135022460 -0.002617965 FESARU
## GALAPA 0.077753255 0.581 -0.0146881661 0.013442529 GALAPA
## GALPED 0.100065668 0.350 0.0224825505 0.019236293 GALPED
## GERCAR 0.009393616 0.898 0.0108559078 -0.010491515 GERCAR
## GEUCAN 0.022772173 1.000 -0.0142451724 -0.012361753 GEUCAN
## HELAUT 0.056859246 0.670 -0.0583436417 -0.086514015 HELAUT
## HELMOL 0.444901394 0.002 0.3553255256 -0.034761103 HELMOL
## HORPUS 0.312087172 0.017 -0.1475521122 -0.293480395 HORPUS
## IPOLAC 0.009490792 0.913 -0.0090387556 -0.031331641 IPOLAC
## IVAANN 0.311525948 0.008 0.0455330518 0.039201868 IVAANN
## JUNINT 0.133690558 0.202 0.0136331062 -0.016777789 JUNINT
## KOEMAC 0.370585711 0.009 -0.0487436827 0.170858654 KOEMAC
## LEPVIR 0.325601019 0.016 -0.0572340552 0.003521346 LEPVIR
## LESCAP 0.334806572 0.019 0.1504633256 -0.021946339 LESCAP
## LIAPYC 0.163347246 0.049 0.0179812798 0.015481057 LIAPYC
## LONMAC 0.172667641 0.105 0.0004196733 0.031553343 LONMAC
## LOTCOR 0.348615977 0.013 -0.0022111948 0.443190652 LOTCOR
## MONFIS 0.091609045 0.361 -0.0331070975 0.064020687 MONFIS
## MYOVER 0.173081368 0.121 -0.0378162987 0.046423163 MYOVER
## OXADIL 0.143002494 0.209 -0.0721611624 0.020946220 OXADIL
## PANDIC 0.006977944 0.921 -0.0207573701 0.008207711 PANDIC
## PENDIG 0.193114410 0.108 -0.1192823538 0.025241468 PENDIG
## PHAARU 0.041539879 0.794 -0.0859544746 -0.074776638 PHAARU
## PLALAN 0.001264554 0.986 0.0202676327 0.017207299 PLALAN
## PLAMAJ 0.132069256 0.231 -0.0162510996 0.014872917 PLAMAJ
## PLAVIR 0.091934859 0.479 -0.0141967833 0.012992818 PLAVIR
## POAPRA 0.137217443 0.219 0.0492069968 0.092727229 POAPRA
## PYCTEN 0.083336781 0.521 -0.0189514955 -0.016445831 PYCTEN
## RANABO 0.181282020 0.125 0.0487820743 0.041907584 RANABO
## RATPIN 0.671445336 0.001 0.5240609418 -0.268029972 RATPIN
## RUDHIR 0.049085637 0.584 -0.0093949298 -0.096825923 RUDHIR
## RUMCRI 0.068456245 0.460 -0.0807045564 0.043057146 RUMCRI
## SETFAB 0.118809057 0.255 -0.0327750423 -0.097865804 SETFAB
## SETPUM 0.545215500 0.001 -0.4458985710 0.039650195 SETPUM
## SIDSPI 0.059385684 0.677 -0.0161814285 -0.014042007 SIDSPI
## SOLCAR 0.007758942 0.916 -0.0193228840 -0.027701254 SOLCAR
## SORNUT 0.789440722 0.001 0.7621760713 0.011031048 SORNUT
## SPHOBT 0.171127596 0.118 0.0372025018 0.032029647 SPHOBT
## SPOHET 0.163347246 0.049 0.0440449605 0.037920690 SPOHET
## SYMPIL 0.105918848 0.330 -0.0346702003 0.027790858 SYMPIL
## TAROFF 0.111769896 0.291 -0.0852401230 -0.004968236 TAROFF
## TRIPRA 0.163347246 0.049 0.0672797885 0.057924811 TRIPRA
## TRIREP 0.412564825 0.003 -0.1045058986 0.087732099 TRIREP
## VERARV 0.005195253 0.958 0.0221452014 0.016543448 VERARV
## VERURT 0.133690558 0.202 0.0136331062 -0.016777789 VERURT
## VIOPED 0.211880088 0.036 -0.0363283830 0.033247537 VIOPED## r p axis1 axis2 labels
## ACHMIL 0.5570346 0.001 -0.07111052 0.59252083 ACHMIL
## CORTRI 0.6359363 0.001 0.48207605 0.11241461 CORTRI
## HELMOL 0.4449014 0.002 0.35532553 -0.03476110 HELMOL
## IVAANN 0.3115259 0.008 0.04553305 0.03920187 IVAANN
## KOEMAC 0.3705857 0.009 -0.04874368 0.17085865 KOEMAC
## RATPIN 0.6714453 0.001 0.52406094 -0.26802997 RATPIN
## SETPUM 0.5452155 0.001 -0.44589857 0.03965020 SETPUM
## SORNUT 0.7894407 0.001 0.76217607 0.01103105 SORNUT
## TRIREP 0.4125648 0.003 -0.10450590 0.08773210 TRIREP
2023 NMDS
adonis2(inner.hel.2023 ~ Treatment, data = env.data.2023, dist="bray")## Permutation test for adonis under reduced model
## Terms added sequentially (first to last)
## Permutation: free
## Number of permutations: 999
##
## adonis2(formula = inner.hel.2023 ~ Treatment, data = env.data.2023, dist = "bray")
## Df SumOfSqs R2 F Pr(>F)
## Treatment 3 1.2223 0.32383 3.1928 0.001 ***
## Residual 20 2.5522 0.67617
## Total 23 3.7745 1.00000
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1### Pair-wise comparisons
library(pairwiseAdonis)
pairwise.adonis2(inner.hel.2023 ~ Treatment, data = env.data.2023, p.adjust = "bonferroni")## $parent_call
## [1] "inner.hel.2023 ~ Treatment , strata = Null , permutations 999"
##
## $LE_vs_LL
## Df SumOfSqs R2 F Pr(>F)
## Treatment 1 0.55758 0.3416 5.1883 0.004 **
## Residual 10 1.07469 0.6584
## Total 11 1.63228 1.0000
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## $LE_vs_NAT
## Df SumOfSqs R2 F Pr(>F)
## Treatment 1 0.22834 0.13662 1.5823 0.106
## Residual 10 1.44310 0.86338
## Total 11 1.67144 1.00000
##
## $LE_vs_SIM
## Df SumOfSqs R2 F Pr(>F)
## Treatment 1 0.44745 0.25705 3.4599 0.003 **
## Residual 10 1.29325 0.74295
## Total 11 1.74070 1.00000
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## $LL_vs_NAT
## Df SumOfSqs R2 F Pr(>F)
## Treatment 1 0.34918 0.21714 2.7736 0.005 **
## Residual 10 1.25892 0.78286
## Total 11 1.60810 1.00000
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## $LL_vs_SIM
## Df SumOfSqs R2 F Pr(>F)
## Treatment 1 0.2803 0.20175 2.5274 0.004 **
## Residual 10 1.1091 0.79825
## Total 11 1.3894 1.00000
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## $NAT_vs_SIM
## Df SumOfSqs R2 F Pr(>F)
## Treatment 1 0.58171 0.28249 3.9372 0.006 **
## Residual 10 1.47748 0.71751
## Total 11 2.05919 1.00000
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## attr(,"class")
## [1] "pwadstrata" "list"Plot Panel
